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ABSTRACT 

Taylor series integration is implemented in NASA Glenn’s Spacecraft N-body Analysis 
Program, and compared head-to-head with the code’s existing 8 th - order Runge-Kutta Fehlberg 
time integration scheme. This paper focuses on trajectory problems that include oblateness 
and/or variable atmospheric density. Taylor series is shown to be significantly faster and more 
accurate for oblateness problems up through a 4x4 field, with speedups ranging from a factor of 
2 to 13. For problems with variable atmospheric density, speedups average 24 for atmospheric 
density alone, and average 1.6 to 8.2 when density and oblateness are combined. 

INTRODUCTION 

We recently demonstrated the significant advantages of Taylor series integration in calculating 
spacecraft trajectories, i.e., order-of-magnitude speedup and simultaneous improvement in accu- 
racy over a conventional 8 th -order Runge-Kutta method. 1 Force models included central body 
gravitation (including J 2 ), other body gravitation, thrust, constant atmospheric drag and solar rad- 
iation pressure with constant illumination. 

In this paper we consider gravitational harmonics beyond J 2 and drag due to variable atmos- 
pheric density. We show that Taylor series (TS) offers significant computational advantages in 
both cases, but also has limitations. 

An alternative TS approach for handling oblateness and variable atmospheric density has been 
presented by Montenbruck. 2 The key differences here are the present implementation applies to 
spacecraft virtually anywhere in the solar system, can be used interchangeably with another inte- 
gration method and employs the 1976 Standard Atmosphere 3 . 

In the next section we summarize the TS methodology for handling oblateness effects and at- 
mospheric density, and discuss its implementation in NASA Glenn’s Spacecraft N-body Analysis 
Program (SNAP) 4 . Following that we present several example problems and discuss numerical 
results. Head-to-head comparisons are made with SNAP’s 8 th -order Runge-Kutta Fehlberg 
(RKF) scheme. 5 
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TAYLOR SERIES FORMULATION 


The governing equations of spacecraft motion are 

x' - v 

v' = a(x 1 ,x 2 ,x 3 ,x 4 ,x s ,x 6 ,x 7 ,t) (1) 

x 7 ' = — m(t ) 

where ( x 1 ,x 2 ,x 3 ) = x is the spacecraft position in Cartesian coordinates relative to an inertial 
frame centered at the central body, (x 4 , x 5 , x 6 ) = v is the velocity, x 7 is the spacecraft mass, 
m is the mass flow rate and a is the acceleration, where 

Cl Cl c b + Cl 0 b T Ctth T Ctd T GE sr p T ^obc T cl o4)0 (2) 

a is the sum of accelerations from the central body point mass, other body point masses, thrust, 
atmospheric drag, solar radiation pressure, oblateness effects of the central body, and oblateness 
effects of other bodies. 

As described previously 1 , the basic approach is to expand each component of the spacecraft 
state vector X = (x 1 ,x 2 ,x 3 ,x 4 ,x 5 ,x 6 ,x 7 ) in a local Taylor series which is then used to advance the 
solution to the next time level. Derivatives are obtained by reformulating Equations (1) into a 
desired canonical form and then making use of highly efficient differentiation arithmetic. A se- 
ries expansion is made at each time level i, 

x,(‘)=£x,m(‘-o‘ , o) 

k = 0 


for n = 1,2, . .., N, where N is the number of equations in the reformulated system, and the Taylor 
coefficients X n (k) are evaluated through simple recursive formulas. 

Oblateness Effects of the Central Body 

The gravitational potential due to a non-spherical shaped central body is expressed in the form 
of a truncated series whose coefficients are obtained from experimental observation. Non- 
spherical gravity accelerations are readily obtained as the gradient of the potential function. 
While simple in concept, the mathematical expressions are quite involved and will not be pre- 
sented here. See Baker 6 (pp. 144-150), Battin 7 (pp. 401-408) or Bate et al 8 (pp. 419-423). 


Table 1. Number of Auxiliary Variables per Oblateness Case 


n 

m 

0 

2 

3 

4 

0 

2 

33 

52 

59 

1 


59 

74 

90 

2 


68 

90 

113 

3 



99 

129 

4 




138 
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Due to the complexity of the oblateness equations, the number of series terms employed in the 
gravitational potential and the need to transform from earth-fixed coordinates to the equator and 
equinox-of-date framework, the number of auxiliary variables required to handle oblateness ef- 
fects is significantly higher than that of other acceleration terms. Table 1 summarizes the number 
of auxiliary variables required up through a 4 x 4 field. Each auxiliary variable corresponds to 
an additional equation in the differential system. The number of auxiliary variables required for 
the central body point mass is shown as the n=0, m=0 case, where (n,m) = (degree, order) of the 
spherical harmonics. Note that the central body point mass requires only two additional equa- 
tions, whereas J 2 alone requires 33 additional equations. The (4,4) case requires 138 additional 
equations. Clearly, the increasing number of equations will make the TS method inefficient as 
oblateness resolution is increased beyond a certain point. 


Variable Atmospheric Density 


The drag acceleration term a d is 



a d,l = Cl [(*4 + C 2 W ) 2 + (*5 - C 2*l ) 2 + A . 

2 (x 4 + c 2 x 2 )/ x 7 

(4) 

Ct d 2 = Cj [fa + C 2 X 2 ) 2 + (x 5 - C 2 X x \ + xj 

2 (x 5 - C 2 Xj)/ x 7 

(5) 

a d p = Cj [(x 4 + c 2 x 2 ) 2 + (x 5 - c 2 x t y + x 6 2 ] 2 x 6 / x 7 

(6) 


where c l = — — C D ■ A ■ p , A is the spacecraft area, p is the atmospheric density, C D is the 


drag coefficient and c 2 = rotation rate of Earth. For spacecraft orbits where p varies, it is ne- 
cessary to calculate the time derivatives of c 2 . The first derivative is just 


dc 2 1 dp 

~dt = ~ 2 C ° A dt 


(7) 


where 


dp dp dp dx 2 dp dx 2 dp dx 3 

dt dt dx 2 dt dx 2 dt dx 3 dt 


( 8 ) 


or 


dp 

dt 


dp 

= -+Vp.„ 


(9) 


The derivatives of p must be obtained from the atmospheric density model. We show how this 
can be done using the 1976 Standard Atmosphere. 

Let the spacecraft altitude h be approximated by 

h= \x\- r e (10) 

where r e is the radius of the earth. Then 


dp dp dh dp x ± 

dx ± dh dx 2 dh |x| 


( 11 ) 
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( 12 ) 


dp dp dh dp x 2 

dx 2 dh dx 2 dh |x| 


so that 


dp dp dh dp x 3 

dx 3 dh dx 3 dh |x| 


dp 1 
Vp = dh \x\ 


x 


Since ^ = 0 for the 1976 Standard Atmosphere, Equation (9) becomes 


dp dp 1 

dt dh |x| 


( 13 ) 


( 14 ) 


( 15 ) 


Equation (15) can be evaluated so long as jj- can be determined for arbitrary h. Here ^ can be 
approximated using finite differences. For h bracketed by h t and h i+1 in the Standard Atmos- 


phere discretization, 
least first order. 


— — approximates ^ to at least first order, so 

h i+1 - hi dh 


dp 

dt 


is also known to at 


Equation (15) shows that ^ is inversely proportional to the magnitude of the position vector. 

In addition, for a purely circular orbit, x ■ v = 0, so that ^ is zero. This is consistent with and 

follows from the fact that density is a function of altitude only in the 1976 Standard Atmosphere. 
Of course, there are no perfectly circular spacecraft orbits, but in general x ■ v will be quite 

small for orbits with small eccentricity. The result is that ^ will be small and higher time deriva- 
tives of density may be negligible in many instances. 


APPLICATION TO NEAR EARTH TRAJECTORIES 

We consider the following problems. (1) Earth-orbiting satellite with oblateness effects and 
zero atmospheric drag; (2) Spacecraft spiraling out of earth’s gravity well in a low thrust trajecto- 
ry with oblateness effects and zero atmospheric drag; (3) Earth-orbiting satellite with atmospheric 
drag and zero oblateness; (4) Earth-orbiting satellite with variable atmospheric drag and oblate- 
ness effects; (5) Earth-orbiting satellite with constant atmospheric drag and oblateness effects. 
For all problems the spacecraft mass is 10,000 kg and the moon and sun are included as other bo- 
dies. The drag coefficient is 0.7 and the spacecraft area is 240 square meters. In Problems (1) 
and (3) - (5), the earth-orbiting satellite is at an inclination of 28.45 degrees with initial altitude 
415 km and the trajectory is propagated for 10 days. For Problem 2 the calculation ends when the 
semimajor axis equals 40,000 km. 

All calculations were run on a Dell Poweredge 6950, which has 8 AMD 2.8 GHz processers 
and 16 GB of RAM. TS calculations used a series with 20 terms for all problems. Both RKF 
and TS used a local error tolerance defined by T{ = Min n [ \x n \ t + x] where x n is the n th varia- 
ble in the state vector/differential system, and x denotes the error tolerance parameter. 
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Problems with Oblateness 


We consider Problems 1 and 2, with spherical harmonics having degree and order (n,m) = 
(2,0), (2,2), (3,3) and (4,4). The tables below show (a) the head-to-heacT ratio of RKF/TS CPU 
time, (b) relative L 2 convergence of each method to its x = l.e-14 solution and (c) difference be- 
tween TS and RKF solutions at x = l.e-14. Relative L 2 convergence is defined by |x — Xf\l |x^| 
where Xf is the final spacecraft position at x = l.e-14 and x is the final spacecraft position at 
larger x, as listed in the tables. Relative L 2 difference between RKF and TS is defined by 
I*ts — *rkfI / I*rkfI- 

Table 2. a CPU ratios, RKF/TS, for Problem 1 


(n,m) 

X 

(2,0) 

(2,2) 

(3,3) 

(4,4) 

l.e-10 

8.61 

3.07 

2.29 

1.61 

l.e-11 

9.87 

3.58 

2.57 

1.89 

l.e-12 

12.94 

4.12 

3.10 

2.22 

l.e-13 

15.14 

4.88 

3.59 

2.68 

l.e-14 

17.41 

5.83 

4.31 

3.17 

Average: 

12.80 

4.29 

3.17 

2.31 


Table 2.b Relative L 2 convergence levels for Problem 1 
RKF TS 


(n,m) 

X 

(2,0) 

(2,2) 

(3,3) 

(4,4) 

(2,0) 

(2,2) 

(3,3) 

(4,4) 

l.e-10 

0.42e-06 

0.42e-06 

0.42e-06 

0.42e-06 

0.13e-07 

0.36e-10 

0.13e-09 

0.24e-10 

l.e-11 

0.34e-07 

0.34e-07 

0.34e-07 

0.34e-07 

0.56e-08 

0.25e-10 

0.71e-10 

0.54e-10 

l.e-12 

0.27e-08 

0.27e-08 

0.26e-08 

0.27e-08 

0.31e-09 

0.79e-10 

0.10e-09 

0.94e-10 

l.e-13 

0.22e-09 

0.21e-09 

0.18e-09 

0.20e-09 

0.35e-09 

0.14e-09 

0.47e-10 

0.85e-10 


Table 2.c 

Difference between RKF and TS solution for Problem 1 

(n,m) 

Relative L 2 Difference 

Absolute Difference (km) 

(2,0) 

0.148e-09 

0.100e-05 

(2,2) 

0.322e-10 

0.218e-06 

(3,3) 

0.819e-10 

0.554e-06 

(4,4) 

0.146e-09 

0.990e-06 


Table 3. a CPU ratios, RKF/TS, for Problem 2 


(n,m) 

X 

(2,0) 

(2,2) 

(3,3) 

(4,4) 

l.e-10 

9.14 

3.23 

2.36 

1.74 

l.e-11 

10.60 

3.74 

2.76 

2.03 

l.e-12 

12.55 

4.45 

3.27 

2.39 

l.e-13 

15.32 

5.23 

3.86 

2.83 

l.e-14 

17.72 

6.16 

4.53 

3.34 

Average: 

13.07 

4.56 

3.35 

2.47 
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Table 3.b Relative L 2 convergence levels for Problem 2 
RKF TS 


(n,m) 

(2,0) 

(2,2) 

(3,3) 

(4,4) 

(2,0) 

(2,2) 

(3,3) 

(4,4) 

T 

l.e-10 

0.13e-04 

0.13e-04 

0.13e-04 

0.13e-04 

0.32e-05 

0.21e-07 

0.15e-08 

0.10e-08 

l.e-11 

0.11e-05 

0.11e-05 

0.11e-05 

0.1 le-05 

0.87e-07 

0.18e-08 

0.32e-09 

0.86e-09 

l.e-12 

0.83e-07 

0.83e-07 

0.83e-07 

0.83e-07 

0.30e-07 

0.77e-09 

0.25e-09 

0.12e-08 

l.e-13 

0.59e-08 

0.58e-08 

0.63e-08 

0.61e-08 

0.58e-08 

0.13e-09 

0.85e-09 

0.16e-09 


Table 3.c 

Difference between RKF and TS solution for Problem 2 

(n,m) 

Relative L 2 Difference 

Absolute Difference (km) 

(2,0) 

0.526e-07 

0.210e-02 

(2,2) 

0.662e-07 

0.265e-02 

(3,3) 

0.753e-07 

0.301e-02 

(4,4) 

0.769e-07 

0.308e-02 


Average TS speedups range from 2.31 to 12.8 for Problem 1, and from 2.47 to 13.07 for Prob- 
lem 2. Tables 2.b and 3.b show that TS is more converged than RKF in all but one case. CPU 
times versus convergence levels are plotted in Figures 1 and 2. 

The small differences in RKF and TS solutions shown in Tables 2.c and 3.c are due to the fact 
that TS calculates the motion of the sun and moon, whereas RKF obtains their motion from 
ephemeris files. The differences are bigger in Problem 2 because the trajectory is propagated 
over 112 days versus 10 days for Problem 1. 



Figure l.a CPU time vs. valid digits for Problem 1, (n,m) = (2,0), (2,2) 
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- Log (Convergence Level) 

Figure l.b CPU time vs. valid digits for Problem 1, (n,m) = (3,3), (4,4) 



Figure 2. a CPU time vs. valid digits for Problem 2, (n,m) = (2,0), (2,2) 
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Figure 2.b CPU time vs. valid digits for Problem 2, (n,m) = (3,3), (4,4) 


Problems with Atmospheric Drag 


Table 4. Effect of Time Derivative Density Term on Results 

Difference in RKF and TS spacecraft position (km) 


Without Variable With Variable 

Density Model Density Model 


l.e-10 0.0306 0.00058 

l.e-11 0.0292 0.00189 

l.e-12 0.0256 0.00153 

l.e-13 0.0230 0.00148 

l.e-14 0.0205 0.00124 


We present results for an earth-orbiting satellite with atmospheric drag but no oblateness 
(Problem 3). We first consider the effect of the time derivative density term Table 4 shows 
the difference in RKF and TS spacecraft positions after the satellite trajectory has been propagat- 
ed for 10 days, where the TS trajectory is propagated with and without . The results clearly 
show the improvement in accuracy due to inclusion of the time derivative density term. 

Table 5 summarizes RKF and TS performance on Problem 3. Constant density results are in- 
cluded for comparison. The speedups are quite large here due to the fact that TS can handle the 
drag equations with only four auxiliary variables. Figure 3 plots the CPU times versus conver- 
gence levels. 
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Table 5. a CPU ratios, RKF/TS, for Problem 3 


I 

Constant Density 

Variable Density 

l.e-10 

23.38 

17.39 

l.e-11 

28.82 

19.85 

l.e-12 

34.84 

23.30 

l.e-13 

40.77 

28.22 

l.e-14 

47.90 

33.41 

Average: 

35.14 

24.43 



Table 5.b Relative L 2 convergence levels for Problem 3 




RKF 

TS 


T 

Constant Density 

Variable Density 

Constant Density 

Variable Density 

l.e-10 

0.43e-06 

0.44e-06 

0.68e-08 

0.17e-06 

l.e-11 

0.35e-07 

0.35e-07 

0.10e-09 

0.13e-06 

l.e-12 

0.27e-08 

0.28e-08 

0.37e-09 

0.46e-07 

l.e-13 

0.20e-09 

0.20e-09 

0.57e-10 

0.36e-07 


Table 5.e Difference between RKF and TS solution for Problem 3 

Constant Density Variable Density 

Relative L 2 Difference Absolute Difference (km) Relative L 2 Difference Absolute Difference (km) 

0.117e-10 0.794e-07 0.182e-06 0.124e-02 



- Log (Convergence Level) 

Figure 3. CPU time vs. valid digits for Problem 3 
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Problems with Atmospheric Drag and Oblateness 

We consider Problems 4 and 5. RKF and TS performance are summarized below and results 
are plotted in the figures which follow. Average TS speedups range from 1.6 to 8.2 for Problem 
4 and from 2.2 to 11.9 for Problem 5. The largest absolute difference between RKF and TS 
shown in Table 6.c is 81.8 meters. It is also worth noting that RKF and TS both predict a reduc- 
tion in the orbit semimajor axis of 2.709 km (for the (4,4) case). This compares with the drag- 
only case (Problem 3) in which both methods predict a reduction in semimajor axis of 1.906 km. 
Atmospheric drag combined with earth’s oblateness thus leads to a significantly greater loss in 
spacecraft altitude. 

Table 6. a CPU ratios, RKF/TS, for Problem 4 




(n,m) 

T 

(2,0) 

(2,2) 

(3,3) 

(4,4) 





l 

l.e-10 

5.89 

2.12 

1.55 

1.13 





l.e-11 

6.84 

2.47 

1.81 

1.32 





l.e-12 

7.96 

2.88 

2.09 

1.53 





l.e-13 

9.19 

3.41 

2.45 

1.84 





l.e-14 

11.03 

3.98 

2.88 

2.17 





Average: 

8.18 

2.97 

2.16 

1.60 





Table 6.b Relative 

L 2 convergence 

levels for Problem 4 





RKF 




TS 


(n,m) 

X 

(2,0) 

(2,2) (3,3) 

(4,4) 


(2,0) 

(2,2) 

(3,3) 

(4,4) 

l.e-10 

0.42e-06 

0.43e-06 0.43e-06 

0.43e-06 


0.40e-05 

0.20e-04 

0.14e-04 

0.12e-04 

l.e-11 

0.35e-07 

0.36e-07 0.34e-07 

0.37e-07 


0.12e-05 

0.13e-04 

0.86e-05 

0.73e-05 

l.e-12 

0.10e-08 

0.36e-08 0.28e-08 

0.48e-08 


0.94e-05 

0.92e-05 

0.44e-05 

0.44e-05 

l.e-13 

0.92e-09 

0.22e-08 0.73e-09 

0.48e-09 


0.77e-05 

0.36e-05 

0.20e-05 

0.19e-05 


Table 6.c Difference between RKF and TS solution for Problem 4 


(n,m) 

Relative L 2 Difference 

Absolute Differe 

(2,0) 

0.881e-06 

0.597e-02 

(2,2) 

0. 121e-04 

0.818e-01 

(3,3) 

0.856e-05 

0.579e-01 

(4,4) 

0.732e-05 

0.495e-01 


Table 7. a CPU ratios, RKF/TS, for Problem 5 


(n,m) 

i 

(2,0) 

(2,2) 

(3,3) 

(4,4) 

l.e-10 

8.14 

2.98 

2.12 

1.62 

l.e-11 

9.44 

3.43 

2.52 

1.84 

l.e-12 

11.87 

4.03 

2.92 

2.18 

l.e-13 

13.70 

4.72 

3.44 

2.53 

l.e-14 

16.42 

5.71 

4.07 

3.05 

Average: 

11.91 

4.17 

3.02 

2.24 
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Table 7.b Relative L 2 convergence levels for Problem 5 
RKF TS 


(n,m) 

T 

(2,0) 

(2,2) 

(3,3) 

(4,4) 

(2,0) 

(2,2) 

(3,3) 

(4,4) 

l.e-10 

0.42e-06 

0.42e-06 

0.42e-06 

0.42e-06 

0.29e-07 

0.13e-09 

0.69e-10 

0.88e-10 

l.e-11 

0.34e-07 

0.34e-07 

0.34e-07 

0.34e-07 

0.39e-08 

0.45e-10 

0.67e-10 

0.39e-10 

l.e-12 

0.27e-08 

0.27e-08 

0.27e-08 

0.26e-08 

0.34e-09 

0.48e-10 

0.39e-10 

0.75e-10 

l.e-13 

0.20e-09 

0.15e-09 

0.21e-09 

0.21e-09 

0.51e-10 

0.14e-09 

0.71e-10 

0.45e-10 


Table 7.c Difference between RKF and TS solution for Problem 5 


(n,m) 

Relative L 2 Difference 

Absolute Differe 

(2,0) 

0.843e-10 

0.570e-06 

(2,2) 

0.282e-10 

0.191e-06 

(3,3) 

0.700e-10 

0.474e-06 

(4,4) 

0.563e-10 

0.381e-06 


u 

<u 

<u 

E 


3 

Cl 

u 



( 2 , 0 ) 

( 2 , 0 ) 

( 2 , 2 ) 

( 2 , 2 ) 


- Log (Convergence Level) 

Figure 4. a CPU time vs. valid digits for Problem 4, (n,m) = (2,0), (2,2) 
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Figure 4.b CPU time vs. valid digits for Problem 4, (n,m) = (3,3), (4,4) 



Figure 5. a CPU time vs. valid digits for Problem 5, (n,m) = (2,0), (2,2) 
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Figure 5.b CPU time vs. valid digits for Problem 5, (n,m) = (3,3), (4,4) 


CONCLUSION 

Taylor series integration is implemented in a high fidelity trajectory propagation code and ap- 
plied to near earth trajectory problems. Head-to-head comparisons are made with 8 th -order 
Runge-Kutta Fehlberg integration. Results show that TS is faster by a factor of 2-13 for oblate- 
ness problems up through a 4x4 field, while simultaneously improving accuracy. For problems 
including both oblateness and variable atmospheric density, TS is faster by a factor of 1.6 to 8.2. 
We conclude that TS is a highly competitive integration method for realistic trajectory problems 
involving oblateness and atmospheric density. 
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